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We construct new, efficient, and accurate high-order finite differencing operators which satisfy sum¬ 
mation by parts. Since these operators are not uniquely defined, we consider several optimization 
criteria: minimizing the bandwidth, the truncation error on the boundary points, the spectral radius, 
or a combination of these. We examine in detail a set of operators that are up to tenth order accu¬ 
rate in the interior, and we surprisingly find that a combination of these optimizations can improve 
the operators' spectral radius and accuracy by orders of magnitude in certain cases. We also con¬ 
struct high-order dissipation operators that are compatible with these new finite difference operators 
and which are semi-definite with respect to the appropriate summation by parts scalar product. We 
test the stability and accuracy of these new difference and dissipation operators by evolving a three- 
dimensional scalar wave equation on a spherical domain consisting of seven blocks, each discretized 
with a structured grid, and connected through penalty boundary conditions. 


I. INTRODUCTION 

Kreiss and Scherer proposed quite some time ago 
dll a powerful way of constructing linearly stable 
schemes for solving evolution partial differential equa¬ 
tions which admit an energy estimate at the continuum, 
through the use of difference operators satisfying sum¬ 
mation by parts (SBP). These operators essentially make 
it possible, up to boundary terms, to derive estimates 
analogous to the continuum ones. While the latter guar¬ 
antee well posedness, their discrete counterparts guar¬ 
antee numerical stability. The boundary terms left af¬ 
ter SBP can be controlled by, for example, orthogonal 
projections dill, penalty terms | 6], or a combination 
of them ||7] (see 0, 9] for a comparison between these 
methods). Furthermore, SBP difference operators and 
penalty techniques have been rather recently combined 
to construct stable schemes of arbitrary high order for 
multi-block simulations G3E1- These are simulations 
where the domain is broken into different sub-domains 
which are "glued" together through an appropriate in¬ 
terface treatment, in this case penalty terms. This semi- 
structured approach allows for non-trivial geometries 
while at the same time ensuring stability for schemes of 
arbitrary high order using derivatives satisfying SBP. 

Systems which have smooth solutions (that is, with- 
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out shocks), such as the Einstein vacuum equations 
(see, for example, ill ), are ideal for using high order 
methods. Furthermore, in numerical relativity one typ¬ 
ically deals with non-trivial topologies and the need for 
smooth boundaries. Although there are proofs for par¬ 
ticular systems in non-smooth domains, proofs of well 
posedness for the initial-boundary value problem for 
general symmetric hyperbolic systems usually require 
smooth outer boundaries IIT. ImIi^I . 

Multi-block domains are also more efficient than 
single-block ones, as they can be chosen to adapt to 
particular situations. For instance, they can be made 
to mimic spherical coordinates which automatically re¬ 
duce the angular resolution at large radii (this allowed, 
for example, studying late time behavior in a rotating 
black hole background in full three-dimensions, placing 
the outer boundary at very larg e distances with mod¬ 
est computational resources fll6f|'). and one can also re¬ 
duce the radial resolution (e.g. logarithmically). Last, 
the need for non-trivial topologies includes the particu¬ 
lar but very important case of black hole excision, where 
the black hole singularity is removed from the compu¬ 
tational domain. In sum, the penalty multi-block ap¬ 
proach combined with high order SBP operators ap¬ 
pears to be promising for simulating Einstein's equa¬ 
tions. 

In principle, modulo the tedious but straightfor¬ 
ward symbolic manipulation algebra needed to con¬ 
struct high order difference operators satisfying SBP, 
one can systematically generate in this way stable multi¬ 
block schemes of arbitrary high order. However, it turns 
out that high order operators satisfying SBP are highly 
non-unique, the higher their order the higher their non¬ 
uniqueness. There is great variation among the proper¬ 
ties of these operators, and for reasons that we discuss 
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below, much care has to be taken in choosing the sten¬ 
cils if explicit time integration schemes are used. One 
approach could involve choosing operators with mini¬ 
mum bandwidth, as they reduce the number of opera¬ 
tions. Unfortunately, in some cases the resulting oper¬ 
ator leads to an amplification matrix with a very large 
spectral radius (which has already been pointed out in 
iTd l and 11^1 1: when using explicit schemes to integrate 
in time, this translates into a very small Courant limit. 
One can do much better by attempting to minimize the 
spectral radius of the complete operator, rather than its 
bandwidth. This in some cases leads to a Courant limit 
two orders of magnitude larger as compared to the mini¬ 
mum bandwidth case. One can sometimes do even bet¬ 
ter: another feature to take into account is the ampli¬ 
tude of the truncation errors at and close to boundaries. 
As we will show, in some cases one can decrease them 
by orders of magnitude while keeping the spectral radius 
small. What we have just briefly discussed is one of the 
goals of this paper, namely to explicitly construct effi¬ 
cient and accurate high order SBP difference operators, 
and compare the above different criteria that can be used 
in their construction. We consider both diagonal and re¬ 
stricted full (non-diagonal) norm based operators; in the 
first (second) case up to order ten (eight) in the interior. 

In many cases of interest, particularly in non-linear 
ones, one might want to add a small amount of arti¬ 
ficial dissipation to the problem. In order not to spoil 
the available energy estimates, the dissipation operator 
has to be negative semi-definite with respect to the SBP 
scalar product. This is not just a technical detail. As 
we will discuss below, in certain cases of interest the 
use of simple dissipation operators that do not satisfy 
this property (e.g. standard Kreiss-Oliger dissipation in 
the interior and no dissipation near boundaries, a choice 
commonly used in some applications) cannot get rid of 
some instabilities, while better dissipation operators do. 
Even if there are no instabilities, a dissipation operator 
that is non-zero close to boundaries is very useful if one 
wants to smooth out aspects of the solution propaga ting 
through the multi-block interfaces. Mattsson et al. Il8il 
have recently presented a way of constructing dissipa¬ 
tion operators that are indeed negative semi-definite for 
arbitrary SBP scalar products, and which extend all the 
way up to the boundaries. Following this prescription 
one can construct, for any difference operator of arbi¬ 
trary high order satisfying SBP, an associated dissipation 
up to the very boundary points in a systematic way. In this 
paper we do so explicitly, for each of the efficient and 
accurate high order derivatives that we present. This is 
the second goal of our paper. 

The third and final goal is to test these new derivative 
and dissipation operators in three-dimensional multi¬ 
block simulations, making use of the penalty method 
to handle interfaces, as described in llOl Hill . Because 
of the challenge involved in achieving stability for very 
high order schemes in the presence of interfaces, multi¬ 
block domains present an ideal setting for testing the 


new derivative and dissipation operators that we here 
construct. While black hole excision is one of our main 
motivations for using multiple blocks, we will report 
here on simulations on a domain that is useful for sce¬ 
narios that do not involve black hole excision, but still 
need a smooth (e.g. spherical) outer boundary. This 
grid structure should be useful for studies of wave phe¬ 
nomena at large distances from the source, gravitational 
collapse, or Friedrich's conformal approach, where the 
spacetime is compactified, and null infinity is brought 
to a finite computational distance (see e.g. [T9|]). In or¬ 
der to isolate testing numerical stability, accuracy, and 
efficiency of the new high order derivative and dissipa¬ 
tion operators from gauge problems and continuum in¬ 
stabilities typically found in many formulations of the 
Einstein equations, we perform these tests in this paper 
using a simpler three-dimensional system —a massless 
scalar field— and we will report on evolutions of the 
Einstein equations elsewhere. 

This paper is organized as follows. In Section [H] 
we introduce our notation, review shortly the penalty 
method, discuss the relative merits of SBP finite dif¬ 
ferencing operators based on diagonal and on non¬ 
diagonal norms, and summarize the construction of dis¬ 
sipation operators of Mattsson et al. jig]. In Section 
EUl we explain the different strategies that we use in 
constructing the new derivative operators, namely their 
bandwidth, their spectral radius, and their truncation 
error. We introduce our example system of evolution 
equations in Section llVl where we also describe the type 
of three-dimensional multi-block domain that we use to 
test the new difference and dissipation operators. We 
describe our new operators corresponding to diagonal 
and restricted full norms in Section 0 discussing their 
properties and comparing their accuracies in numerical 
tests. Finally, Section fVTl closes with some remarks about 
possible future research directions. 

II. SBP DERIVATIVE AND DISSIPATION OPERATORS 
WITH A HIGH ORDER OF ACCURACY 

A. SBP and penalties 

In this subsection we briefly summarize Section 2 of 
m. We do so essentially to fix our notation, for more 
details see that reference. Consider a computational 
domain [a,b] and a discrete grid consisting of points 
i = 1... n and grid spacing h. A difference operator D is 
said to satisfy SBP on that domain with respect to a pos¬ 
itive definite scalar product E (defined by its coefficients 
Vij) 

n 

(u, v) = h UiVjCij , (1) 

y=i 

if the property 

(u, Dv) + {v, Du) = (uv) \ b 


( 2 ) 
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holds for all grid functions u and v. Similar definitions 
can be introduced for two (and higher) dimensional do¬ 
mains. The scalar product or norm is said to be diagonal 
if 


• Negative A: 

S r = —A + S, Si = 5, with 5 > ^ (9) 


(Tij — Ciifiij / (3) 

that is, if Ujj is diagonal. It is called restricted full if 


a ibi ~ ^hh^hi' 


(4) 


that is, if o~ij is diagonal on the boundary, but may be 
non-diagonal (full) in the interior, ij, £ { I, n } denote 
boundary point indices. 

We now briefly highlight through a simple example 
the main features of the penalty method for multi-block 
evolutions, for more details see EH. We assume that 
the norm is either diagonal or restricted full, since these 
are the cases we actually consider later in this paper. 

The simple example we wish to consider is the advec- 
tion equation. 


d t u = A d x ii 


in the spatial interval (—oo, +oo) with appropriate fall- 
off conditions at infinity, and two grids: a left grid cover¬ 
ing (—oo,0], and a right grid covering [0, +oo). We refer 
to the grid function u on each grid by u l and u r , corre¬ 
sponding to the left and right grids, respectively. The 
problem is discretized using grid spacings It 1 , h 1 on the 
left and right grids —not necessarily equal— and dif¬ 
ference operators D\ D r satisfying SBP with respect to 
scalar products given by the weights crfa r on their indi¬ 
vidual grids. That is, these scalar products are defined 
through 


Jv r - 
hi 1 


0 +0O 

(u l ,v l )=h l £ aijujvj , (u\v r ) = h r Y J Vi i u r i 

i,j=—oo i,j =0 

The semi-discrete equations are written as 


(5) 


3f«i = AD?M i + § 2 r( M o- M o)' (6) 

n cr QQ 

d tU r = A DX + j^(u l 0 -u r 0 ). ( 7 ) 

n ^00 

In the fully non-diagonal case the treatment is slightly 
more complicated, therefore we consider here only the 
diagonal and the restricted full cases. 

One can derive an energy estimate and therefore guar¬ 
antee stability if two conditions are satisfied. One of 
them is A + S r — Sj = 0. The other one imposes an ad¬ 
ditional constraint on the values of S; and S r : 

• Positive A: 

Sj = A + 5, S r = 5, with 6 > — ^ (8) 


• Vanishing A: this can be seen as the limiting case 
of any of the above two. 

For the minimum values of S allowed by the above in¬ 
equalities the energy estimate is the same as for a single 
grid (that is, as if the interface did not exist), while for 
larger values of 5 there is damping in the energy which 
is proportional to the mismatch at the interface. 


B. Diagonal versus non-diagonal norms 

There are several advantages in using one¬ 
dimensional (ID) difference operators satisfying 
SBP with respect to diagonal norms. One of them is 
related to the fact that SBP is guaranteed to hold in 
several dimensions by simply applying the ID operator 
along each direction f|, J !5]. Another advantage is 
related to the following: in order to ensure stability 
through an energy estimate, in many cases one has to 
be able to bound the norm of the commutator between 
the difference operator and the principal part of the 
equations for all resolutions, and this is guaranteed to 
hold in the diagonal case. Finally, the expressions of the 
operators are also somewhat simpler when compared 
to non-diagonal ones. The disadvantage, on the other 
hand, is that the order of accuracy at and close to 
boundaries is half of that in the interior, while in the 
restricted full case the operators lose only one order 
near boundaries mil. 

Though more efficient, schemes based on non¬ 
diagonal norms might have stability problems in the ab¬ 
sence of dissipation. First, it is not guaranteed that SBP 
holds in several dimensions if a difference operator sat¬ 
isfying SBP in ID is applied along each direction. Sec¬ 
ond, there can be problems even in one dimension when 
the system has variable coefficients, since the bounded¬ 
ness of the commutator discussed above is not guaran¬ 
teed to hold either. This is not a feature inherent to fi¬ 
nite difference (FD)-based schemes: the boundedness of 
such commutator is, for example, in general not guar¬ 
anteed either for pseudo-spectral methods in the ab¬ 
sence of filtering, even in periodic domains, when the 
system has variable coefficients ||bij]. Therefore, both in 
the case of pseudo-spectral methods and non-diagonal 
norm based FD schemes, one is in general unable to 
guarantee stability in more than one dimension, or even 
in ID in the variable coefficient case, without filtering or 
dissipation (see also |22]). In the problem at hand, the 
question then is whether one can stabilize the scheme 
through artificial dissipation, without introducing an 
excessive amount of it. Below we will address this ques¬ 
tion in detail, as well as compare diagonal based opera¬ 
tors to their non-diagonal counterparts. 
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In the diagonal case we will consider difference oper¬ 
ators of order two, four, six, eight, and ten in the interior 
(and therefore order one, two, three, four and five, re¬ 
spectively, at and close to boundaries) and denote them 
by D 2 _i,D 4 _ 2 ,D 6 _ 3 , D 8 _ 4 , D 10 _ 5 . In the non-diagonal 
(restricted full) case, we will consider operators of or¬ 
der four, six and eight in the interior (and therefore 
order three, five and seven, respectively, at and close 
to boundaries), and denote them by D 4 _ 3 , Dg_ 5 , Dg_ 7 . 
These operators in general are not unique. For exam¬ 
ple, in the second order in the interior case there is a 
unique operator satisfying SBP, and its norm is diago¬ 
nal, the operator being what we called D 2 _i. With re¬ 
spect to higher order operators, the following holds for 
the diagonal norm based ones: D<j _ 2 is unique, while 
D 6 _ 3 , Dg _4 and D 40-5 comprise a one-, three-, and ten- 
parameter family, respectively. In the restricted full case, 
D 4 _ 3 , D( 3_5 and Dg _7 have three, four and five free pa¬ 
rameters, respectively. 


C. Dissipation operators 

As pointed out in 0 , adding artificial dissipation 
may lead to an unstable scheme unless the dissipation 
operator is compatible with the SBP derivative operator. 
In that reference, the authors present a prescription for 
constructing such operators, which we follow here. In 
short, a compatible dissipation operator, of order 2 p in 
the interior, is constructed as 

A 2p =-ah 2 P (10) 

where a is a positive constant, E is the scalar product 
used in the construction of the SBP operator, and Dp 
is a consistent approximation of dP / dx>’ with minimal 
width 1 . Bp is the so-called boundary operator. The bound¬ 
ary operator is positive semi-definite and its role is to al¬ 
low boundary points to be treated differently from inte¬ 
rior points. Bp cannot be chosen freely, but has to follow 
certain restrictions which we explain below. 

For the diagonal norm operators, choosing Bp to be 
the unit matrix is sufficient to obtain the required pth 
order accuracy near the boundary, which is the same ac¬ 
curacy as the derivative operator. 

In the case of restricted full norm operators, the accu¬ 
racy requirement near the boundary is stricter. The dis¬ 
sipation operator should have order 2 p — 1 at the bound¬ 
ary and order 2 p in the interior, which requires a differ¬ 
ent choice of Bp. We again follow 0 and choose Bp 
to be a diagonal matrix, where the diagonal is the re¬ 
striction onto the grid of a piecewise smooth function. 
The numerical domain is divided into three regions in 


1 "Minimal width" means that the stencil must contain as few points 

as possible. 


each dimension; an interior part and on either side two 
transition regions containing the boundaries. The tran¬ 
sition region has a fixed size that is independent of the 
resolution. Within the transition region the function. Bp, 
increases from 0 (/z p_1 ) at the outer boundary to a con¬ 
stant value 1 at the boundary with the interior region in 
such a way that the derivatives of Bp up to order p 2 
vanish at either ends. In the interior region the function 
has the constant value 1 . 

For the D 4-3 operator. Bp has the value li at the 
boundary and increases linearly to 1 in the transition re¬ 
gion. For the Dg _5 operator, we use a cubic polynomial 
with vanishing derivatives at either end of the transition 
region to increase the value of Bp from h 2 at the bound¬ 
ary to 1 in the interior. For the Dg _7 operator, the bound¬ 
ary values for the transition region are h 3 and 1 , and we 
use a fifth order polynomial to make the first and second 
derivatives vanish at either end of the transition region. 

For the constant a we make the choice a = 2 - 2 P / since 
then the parameter used to specify the strength of the 
dissipation has approximately the same allowed numer¬ 
ical range, independently of the order of the operator. 

Note that in the diagonal case up to order eight in 
the interior, the scalar product E is independent of the 
free parameters, so for a given order the same dissipa¬ 
tion operator is used for all the different operators we 
construct below, while for the higher order diagonal op¬ 
erators and the restricted full norm operators a unique 
dissipation operator has to be constructed for each pa¬ 
rameter choice. 


III. OPTIMIZATION CRITERIA 


We start by fixing some notation. If the accuracy of 
the difference operator D in the interior is 2 p, then there 
are b points at and near the boundaries where the or¬ 
der of D is only q. In the diagonal case one has q = p, 
and in the restricted full case it is q = 2p — 1. We call 
b the boundary width. The difference operator at these 
b points uses (up to) s points to compute the derivative. 
We call s the boundary stencil size. 

When expanding D in a Taylor series one has 


Du \ Xi 


du .„ d^ +1 u 
_ -L r -m _ 

dx 1 (M + 1 

x l 


for i = 1 ... b 


( 11 ) 


where h is the grid spacing and x, = ih. We call c,- the 
error coefficients. 

In what follows, we consider three cases for each fam¬ 
ily of operators of a given order, denoted by: 

• Minimum bandwidth: If there are n free parame¬ 
ters, it is always possible to set n of the derivative 
coefficients to zero, thereby minimizing the com¬ 
putational cost of evaluating the derivatives in the 
boundary region. 
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• Minimum spectral radius: In this case, we calculate 
numerically the eigenvalues of the amplification 
matrix for a test problem, and choose the parame¬ 
ters to minimize the largest eigenvalue. 

• Minimum ABTE: We minimize the average bound¬ 
ary truncation error (ABTE), which we define as 


ABTE := 



( 12 ) 


The test problem that we use to compute the spec¬ 
tral radius of the amplification matrix is the same one 
that was used in 11611 : an advection equation propa¬ 
gating in a periodic domain. Periodicity is enforced 
through an artificial interface boundary via penalties. 2 
We use a penalty parameter 3 = —1/2 (see Section Hit, 
which means that the semi-discrete energy is strictly 
preserved, and that the amplification matrix is anti- 
Hermitian, and therefore the real part of all eigenvalues 
is zero. 3 

Another option na would be to compute the spectral 
radius of the discrete difference operator itself. In this 
case, the spectrum is in general not purely imaginary, 
since the boundary conditions have not been imposed 
yet. In practice we have found, though, that both ap¬ 
proaches lead to similar operators, in the sense that a 
derivative operator with small spectral radius usually 
also leads to amplification matrices for the above test 
problem with small spectral radii as well. 

It is worth pointing out that for the diagonal operator 
case the bandwidth and the ABTE can be globally min¬ 
imized by analytically choosing the parameters, since 
the ABTE is a quadratic function of the parameters and 
therefore has a global minimum. This is not the case for 
the spectral radius. Therefore, when we refer to mini¬ 
mizing the spectral radius, we perform a numerical min¬ 
imization and do not claim that we have actually found 
a global minimum. 


IV. EXAMPLE EVOLUTION SYSTEM AND 
MULTI-BLOCK DOMAIN SETUP 


We test each of the new derivative operators and their 
associated dissipation operators through 3D multi-block 
simulations. In these simulations we solve the scalar 
wave equation 


dtt <p = A <p 


(13) 


in static, curvilinear coordinates. This equation can be 
reformulated as 


df <p 

= II 

(14) 

d f n 

= 7 _1/2 D; (V /2 7 ‘idj'j 

(15) 

df di 

= Di n. 

(16) 


where jy is the Euclidean metric in curvilinear coordi¬ 
nates, 7 = det 7 (/ its determinant, and j l t its inverse. 
Here we use the Einstein summation convention, im¬ 
plicitly summing over repeated indices. 

The advantage of this particular form of the equations 
is that it defines a strictly stable discretization in the di¬ 
agonal case, in the sense that one can show by just using 
SBP (i.e., without needing a discrete Leibniz rule, which 
in general does not hold) that there is a semi-discrete 
energy E which is preserved in time for any difference 
operator D satisfying SBP . 4 This energy is given by 

E = \J ( n2 + f’didjj 7 1/2 dV. (17) 

The geometry of the computational domain in the 
three-dimensional (3D) simulations that we show be¬ 
low is the interior of a sphere. In order to avoid the 
singularities at the origin and poles that spherical co¬ 
ordinates have, we cover the domain with seven blocks: 
surrounding the origin we use a Cartesian, cubic block, 
which is matched (at each face) to a set of six blocks 
which are a deformation of the cubed-sphere coordi¬ 
nates used in [T 6 j]. Figure^ shows an equatorial cut of 
our 3D grid structure, for ll 3 points on each block. The 
blocks touch and have one grid point in common at the 
faces. The grid lines are continuous across interfaces, 
but in general not smooth. 

Each block uses coordinates a, b, c. For the inner block 
these are the standard Cartesian ones: a = x,b = y,c = 
z, while for the six outer blocks they are defined as fol¬ 
lows. The "radial" coordinate c £ [—1,1] is defined by 
inverting the relationship 

r = ^[ r o(l-c) + r 1 (l + c)] (18) 

(where r = x 2 +y 2 + z 2 ). The "angular" coordinates, 
a, b G [— 1 , 1 ], are in turn defined through 

• Neighborhood of positive x axis: 
x = r/F, y = rb/F, z = ar/F 

• Neighborhood of positive y axis: 
x=—br/F, y = r/F, z = ar/F 


2 Truly periodic domains (that is, without an interface) do not require 
boundary derivative operators, and therefore do not constitute a 
useful test here. 

3 As a side remark: we actually compute the eigenvalues of the am¬ 

plification matrix multiplied by the grid spacing h. 


4 That is, the energy is preserved modulo boundary conditions, which 
can inject or remove energy from the system. For example, if maxi¬ 
mally dissipative boundary conditions are used, the energy actually 
decreases as a function of time. 
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Figure 1: An equatorial cut of the 3D multi-block structure 
used in the simulations of this paper. 


• Neighborhood of negative x axis: 
x=—r/F, y = —rb/F, z = ar/F 

• Neighborhood of negative y axis: 
x = br/F, y=—r/F, z = ar/F 

• Neighborhood of positive z axis: 

x = —ar/F, y = rb/F, z = r/F 

• Neighborhood of negative z axis: 
x = ar/F, y = rb/F, z =—r/F 

where r is written in terms of c through M St . and 

p -= / (r 1 -r) + (r-r 0 )E \ 1/2 
V n-r 0 ) 

with E = 1 + a 2 + b 2 . The surface c = 1 corresponds to 
the spherical outer boundary (of radius r = i '\), while 
c = —1 corresponds to the cubic interface boundary 
matched to a cube of length 2r () on each direction. The 
grid structure of FigureQJcorresponds to tq = 1, = 3. 

To get an accurate measure of numerical errors, we 
evolve initial data for which there exists a simple ana¬ 
lytic solution of (1141 for all times, against which we can 
compare the numerical results. We choose initial data 

(p(t = 0) = Acos(2nk-x) (20) 

II(f = 0) = —27rA|k|sin(27rfc-x) ( 21 ) 

dj(t = 0) = -2nAkism(2nk-x). (22) 

The analytic solution for this setup is a plane wave with 
constant amplitude A traveling through the grid in the 
direction of the vector k. In all the simulations that we 


present below we use A = 1.0 and k = (0.2,0.2,0.2), 
i.e., the wave is traveling in the direction of the main 
diagonal. Consistent and stable outer boundary condi¬ 
tions are imposed through penalty terms by penalizing 
the incoming characteristic modes with the difference to 
the exact solution, as introduced in | 6 j]. As mentioned 
above, we set r o = 1 and r\ = 3 in our block system, 
placing the outer boundary at R = 3. 

We use resolutions from h = Aa = Ab = Ac = 0.1 (21 3 
grid points per block) up to h = 0.0125 (161 3 grid points 
per block). Our highest resolution corresponds to about 
400 grid points per wave length. (This figure depends 
on which part of which block one looks at, since the 
wave does not propagate everywhere along grid lines.) 
We use that many grid points per wave length —or, 
equivalently put, we use such a large wave length— be¬ 
cause we are interested in high accuracy. Decreasing the 
wave length is resolution-wise equivalent to using fewer 
grid points per block, which is a case we study with our 
coarse resolution. It would be interesting to study the 
effect of very short length features onto the stability of 
the system, i.e., to use a wave length that cannot be well 
resolved any more. Discrete stability guarantees in this 
case that the evolution remains stable, and we assume 
that a suitable amount of artificial dissipation can help 
in the non-linear case when there is no known energy es¬ 
timate. The size of the time step is chosen to be propor¬ 
tional to the minimal grid spacing in local coordinates 
At = A min( Aa, Ab, Ac) with the Courant factor A. Un¬ 
less otherwise stated, we use A = 0.25. For the penalty 
terms, we used <5 = 0 everywhere. For all the runs with 
dissipation the strength was chosen to be e = 0.4 (see 
Appendix |^2J. For the dissipation operators based on 
a non-diagonal norm, we choose the transition region to 
be 30% of the domain size. Note that except for the Dg _5 
operator (where dissipation is essential in order to sta¬ 
bilize the operator) the differences between results with 
and without dissipation are so small, that we only show 
the results obtained without dissipation. In all cases we 
calculate the error in the numerical solution for <p(x,t) 
with respect to the exact analytical solution. 

We have implemented our code in the Cactus frame¬ 
work 123. 24] using the Carpet infrastructure 1255, 2(5]. 


V. OPERATORS 

A. Operators based on diagonal norms 

1. Operator properties 

We consider first operators that are based on a diag¬ 
onal norm, since this is the easier case. We examine 
here the operators D 6 _ 3 , D 8 _ 4 , and D 10 _ 5 . These opera¬ 
tors have 6, 8, and 11 boundary points, respectively, and 
their maximum stencil sizes are 9,12, and 16 points, re- 
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Table I: Properties of the D 2 _ 4 operator. 


Operator 

Unique 

Spectral radius 

1.414 

ABTE 

0.25 

Cl 

0.5 


Table II: Properties of the D 4 _ 2 operator. 


Operator 

Unique 

Spectral radius 

1.936 

ABTE 

0.2276 

Cl 

-0.4215 

C2 

0.1666 

C3 

-0.0193 

C 4 

-0.037 


spectively . 5 The operators D 2 _i and D 4 _ 2 are also based 
on a diagonal norm. They are unique and have been ex¬ 
amined in 111. For completeness we list their properties 
here as well. 

The D 2 _ i operator formally has two boundary points 
and a maximal stencil size of three points. However, the 
stencil for the second boundary point is the same as the 
interior centered stencil, so in practice it has only one 
boundary point with a stencil size of two points. The 
spectral radius and error coefficient are listed in table 
ID The D 4 _ 2 operator has four boundary points and a 
maximal stencil size of six points. We list its properties 
in table Ini 

The family of Dg _3 operators has one free parame¬ 
ter. The resulting norm is positive definite, and is in¬ 
dependent of this parameter, hence the parameter can 
be freely chosen. For this operator there are very small 
numerical differences in the spectral radius and in the 
truncation error coefficients between the three different 
cases where the bandwidth, the spectral radius and the 


Table III: Properties of the diagonal norm Dg _3 operators. 


Operator 

Minimum 

bandwidth 

Minimum 
spectral radius 

Minimum 

ABTE 

Spectral radius 

2.1287 

2.1077 

2.1082 

ABTE 

0.2716 

0.2563 

0.2558 

Cl 

0.5008 

0.5436 

0.5374 

C2 

-0.1854 

-0.2340 

-0.2270 

C3 

-0.2144 

0.0012 

-0.0300 

c 4 

0.3067 

0.1977 

0.2135 

C 5 

-0.1288 

-0.0546 

-0.0654 

C6 

-0.0286 

-0.0419 

-0.0400 


5 We expected the Djo -5 operator to have 10 boundary points and a 
maximum stencil size of 15 points, but this did not result in a posi¬ 
tive definite norm. 


Table IV: Properties of the diagonal norm Dg_ 4 operators. 


Operator 

Minimum 

bandwidth 

Minimum 
spectral radius 

Minimum 

ABTE 

Spectral radius 

16.0376 

2.229 

2.231 

ABTE 

1.2241 

0.3993 

0.3474 

Cl 

-0.5878 

-0.8277 

-0.8086 

C 2 

0.1068 

0.3682 

0.3439 

C 3 

3.1427 

-0.3819 

0.0228 

c 4 

-0.7918 

-0.2186 

-0.3086 

C 5 

0.9886 

-0.3412 

0.0225 

C 6 

0.3304 

0.3619 

0.2970 

C7 

-0.1995 

-0.1097 

-0.0823 

C 8 

-0.0211 

-0.0465 

-0.0497 


average boundary truncation error are respectively min¬ 
imized, as can be seen in table Hill 

Based on those small differences one would expect 
that there should not be much of a difference in terms of 
accuracy among these three different cases in practical 
simulations. However, it turns out that the minimum 
bandwidth operator in practice leads to very different 
solution errors compared to the minimum ABTE oper¬ 
ator, and that the latter is to be preferred. We did not 
implement the minimum spectral radius operator, since 
the difference in spectral radius is minimal. 

The family of D §_ 4 operators has three free parame¬ 
ters; as in the previous case, the norm is positive definite 
and independent of the parameters, which can therefore 
be freely chosen. We again investigate the properties of 
the operators obtained by minimizing the bandwidth, 
the spectral radius, and the ABTE. Interestingly, we find 
out that there is a one parameter family of operators 
that minimizes the ABTE. Therefore, in the minimum 
ABTE case we make use of this freedom and also de¬ 
crease the spectral radius as much as possible. The re¬ 
sults are shown in tablellVI 

In this case, the minimum bandwidth operator is 
quite unacceptable due to its large spectral radius. (For 
this reason we did not even implement it.) The error co¬ 
efficients are also quite large compared to the two other 
cases. The differences in error coefficients between the 
minimum spectral radius operator and the minimum 
ABTE operator might appear quite small, but as demon¬ 
strated in figure [5|below, there is almost of factor of two 
in the magnitude of the error in our numerical tests (see 
below). 

The family of D 4 q _5 operators turns out to be differ¬ 
ent from the lower order cases. The conditions the SBP 
property impose on the norm do not yield a positive def¬ 
inite solution with a boundary width of 10 points. When 
using a boundary width of 11 points instead, there is 
a free parameter in the norm, which for a very narrow 
range of values does allow it to be positive definite. We 
choose this parameter to be approximately in the center 
of the allowed range. With the larger boundary width, 
there are 10 free parameters in the difference operator. 
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Table V: Properties of the diagonal norm D 10 _5 operators. 


Operator 

Minimum 

bandwidth 

Minimum 

ABTE 

Spectral radius 

995.9 

2.240 

ABTE 

20.534 

0.7661 

Cl 

0.7270 

2.0379 

C2 

0.3034 

-0.8545 

C3 

-10.5442 

-0.0898 

c 4 

0.2690 

0.7250 

C5 

4.8373 

0.2889 

C6 

-62.3907 

-0.0154 

C7 

-1.3649 

-0.6429 

C8 

24.6628 

0.0293 

C 9 

0.1045 

0.7073 

C 10 

-0.4197 

-0.1411 

C 11 

-0.0243 

-0.1469 


Fixing these to give a minimal bandwidth operator re¬ 
sults in an operator with a large ABTE (20.534) and very 
large spectral radius (995.9) that is not of practical use. 
Minimizing the spectral radius in the full 10 dimen¬ 
sional parameter space turned out to be very difficult 
because the largest imaginary eigenvalue does not vary 
smoothly with the parameters. Instead we attempted 
to minimize the average magnitude of all the eigenval¬ 
ues, however the resulting operator turned out to have a 
slightly larger spectral radius than the minimum ABTE 
operator considered next and were therefore not imple¬ 
mented and tested. Minimizing the ABTE instead fixes 
four of the ten parameters and the remaining six can 
then be used to minimize the spectral radius. This re¬ 
sults in an operator with a ABTE of 0.7661 and spectral 
radius of 2.240. When used in practice, it turns out that 
this operator at moderate resolutions has rather large er¬ 
rors compared to the corresponding Dg _4 case. The er¬ 
rors can be reduced by a factor of about 3 by adding 
artificial dissipation, indicating that the errors, though 
not growing in time, are dominated by high frequency 
noise from the boundary derivative operators. Only at 
the highest resolution considered in the tests of this pa¬ 
per is there an advantage in using the Djq-s operators. 
We therefore do not pursue them further in this paper, 
though we might consider their use in other applica¬ 
tions if we need higher resolutions. 

For completeness we list the properties of the mini¬ 
mum bandwidth and ABTE Diq -5 operators in table M 


2. Numerical tests 

In the diagonal case we do not need to add dissipation 
to the equations which we solve here, since our semi¬ 
discrete discretization is strictly stable, as discussed in 
Section GYI This means that the errors cannot grow as 
a function of time at a fixed resolution (i.e., at the semi¬ 
discrete level). However, following 11811 we have con¬ 


structed corresponding dissipation operators for these 
derivatives for future use in non-linear problems |27]. 

a. The operator D 6 _ 3 . Figure |2] shows the results of 
convergence tests in the Loo norm for the D (l _3 case, for 
both the minimum bandwidth and the minimum ABTE 
operators. We define the convergence exponent m as 


m 



(23) 


where Li and E 2 are solution errors and h\ and hj the 
corresponding resolutions. In the minimum bandwidth 
case the convergence exponent gets close to three as res¬ 
olution is increased, i.e., the order is being dominated by 
boundary (outer, interface, or both) effects. On the other 
hand, one can see from the figure that in the minimum 
ABTE case we do get a global convergence exponent 
that gets quite close to four when resolution is increased. 
Figure [3] shows an accuracy comparison between these 
two operators, for the coarsest and highest resolutions 
used in the previous plots, displaying the errors with re¬ 
spect to the exact solution in the Loo norm. The improve¬ 
ment is quite impressive: for the highest resolution that 
we used the error with the minimum ABTE operator is 
around tzvo orders of magnitude smaller. 

We conjecture that this is caused by larger truncation 
errors for this operator near boundaries. This is un¬ 
fortunately not immediately evident when looking at 
the boundary error coefficients c, or the ABTE. How¬ 
ever, three of the inner four error coefficients have ab¬ 
solute values less than 0.1 for the optimized operator, 
whereas this is the case for only 1 error coefficient for 
the standard operator. The fact that the accuracy is also 
higher with the optimized operator also points to the 
fact that the standard operator introduces somehow a 
much larger error. 

As a summary, the minimum ABTE operator is the 
preferred choice in this case, the minimum bandwidth 
Dg _3 operator does not have a large spectral radius in 
comparison, but it does have much larger truncation er¬ 
ror coefficients. 

b. The operator Dg_ 4 . Figure 0] shows the results of 
similar convergence tests, also in the Loo norm, for the 
D 8 _4 operators. As discussed in the previous Section, 
in this case the minimum bandwidth operator has both 
very large spectral radius and truncation error coeffi¬ 
cients, so large that it is actually not worthwhile pre¬ 
senting here details of simulations using it (they actu¬ 
ally crash unless a very small Courant factor is used, 
as expected). In references m and m two different 
sets of parameters were found, both of which reduced 
the spectral radius by around one order of magnitude, 
when compared to the minimum bandwidth one. Here 
we concentrate on comparing an operator constructed 
in a similar way (with slightly smaller spectral radius 
than the ones of 1 116111711 ) —that is, minimizing the spec¬ 
tral radius— with the minimum ABTE operator. We see 
that in both cases we find a global convergence exponent 
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Minimum bandwidth 6-3 operators 



Minimum ABTE 6-3 operators 



Minimum spectral radius 8-4 operators 



Minimum ABTE 8-4 operators 



Time 


Figure 2: Convergence exponents for the minimum bandwidth Figure 4: Convergence exponents for the minimum spectral 
(top) and the minimum ABTE (bottom) Dg _3 operators. radius (top) and the minimum ABTE (bottom) Dg _ 4 operators. 
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Figure 3: Accuracy comparison between the Dg _3 operators of 
the previous figure. For the highest resolution that we used, 
the errors with the minimum average boundary truncation er¬ 
ror operator are around two orders of magnitude smaller than 
those obtained with the minimum bandwidth one. 


close to five. Figure [5] shows at fixed resolution (with 
the highest resolution that we used for the convergence 
tests) a comparison between these two operators, by dis¬ 
playing the errors with respect to the exact solution, in 
the Loo norm. There is an improvement of a factor of 
two in the minimum ABTE case (as mentioned, the dif¬ 
ferences with the minimum bandwidth case are much 
larger). Notice also that even though not at round-off 
level, the errors in our simulations are quite small, of 
the order of 10 7 in the Loo norm (in the L 2 norm they 
are almost an order of magnitude smaller). 


B. Operators based on a restricted full norm 

1. Operator properties 

Let us now consider operators that are based on a re¬ 
stricted full norm (see Section ITl At . In this case the norm 
always depends on the free parameters, and it is not nec¬ 
essarily positive definite for all values of them. There¬ 
fore these free parameters are subject to the constraint 
of defining a positive definite norm. 

We examine here the operators D 4 _ 3 , Dg_ 5 , and Dg_ 7 . 


6-3 operators 
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8-4 operators 



Figure 5: Comparison of the accuracy of the two D 8 _ 4 operator 
types shown in the previous figure in the Lnorm. Although 
both operators have quite similar error coefficients, there is up 
to a factor of two difference in the errors seen in the actual 
runs. 


Table VI: Properties of the restricted full norm D 4 _ 3 operators. 


Operator 

Minimum 

bandwidth 

Minimum 
spectral radius 

Minimum 

ABTE 

Spectral radius 

2.428 

1.322 

2.758 

ABTE 

0.1281 

0.5824 

0.0230 

Cl 

0.2500 

1.2359 

- 0.0100 

C2 

-0.1333 

-0.2378 

-0.0409 

C3 

0.0201 

-0.3166 

0.0251 

c 4 

0.0366 

0.1053 

-0.0086 

C5 

-0.0065 

0.0279 

-0.0129 


These operators have five, seven and nine boundary 
points, respectively, and their maximum stencil size are 
seven, ten and thirteen points, respectively. 

The family of D 4 _ 3 operators has three independent 
parameters, and as mentioned above, they have to be 
chosen so that the corresponding norm is positive defi¬ 
nite. 

We constructed operators by minimizing their band¬ 
width, their spectral radius, and their average boundary 
truncation error. In minimizing the bandwidth there is 
some arbitrariness in the choice as to which coefficients 
in the stencils are set to zero. Here we follow 120] and set 
the coefficients 9 / 1 , 5 , q 2 , 7 / and c /3 7 to zero, where the first 
index labels the stencil starting with 1 from the bound¬ 
ary and the second index labels the point in the stencil. 
This results in two solutions, but only one of them cor¬ 
responds to a positive definite norm. 

A comparison of the spectral radius, ABTE, and error 
coefficients is listed in table lVIl 

We find in our simulations (see below) that the min¬ 
imum bandwidth and the minimum ABTE D 4 _ 3 oper¬ 
ators have very similar properties. The latter, however, 
has slightly smaller errors, enough to offset the slightly 
smaller time steps required for stability. The operator 


Table VII: Properties of the restricted full norm Dg _5 operators. 


Operator 

Minimum 

bandwidth 

Minimum 
spectral radius 

Minimum 

ABTE 

Spectral radius 

2.940 

1.458 

3.194 

ABTE 

0.0986 

0.5380 

0.0648 

Cl 

0.1667 

1.3692 

-0.0154 

C2 

-0.1558 

-0.2682 

-0.0507 

C 3 

0.0672 

-0.2118 

0.1336 

c 4 

0.0953 

0.0097 

0.0532 

C 5 

-0.0433 

0.0702 

-0.0733 

c 6 

0.0141 

0.1434 

0.0187 

C7 

-0.0163 

-0.0972 

-0.0123 


with a minimum spectral radius unfortunately has very 
large errors; in fact, we have not been able to stabi¬ 
lize the system with any amount of artificial dissipa¬ 
tion when we used this operator for equations with non¬ 
constant coefficients in 3D. 

The family of Dg_ 5 operators has four independent 
parameters that again have to be chosen so that the 
norm is positive definite. For the minimal bandwidth 
operator we choose to zero the coefficients 971 , 7 , 1/2,9, 9/2,10 
and q 3W . The resulting equations cannot be solved an¬ 
alytically, but numerically we find eight solutions, of 
which four are complex. From the remaining real so¬ 
lutions only one of them results in a positive definite 
norm. A comparison between the properties of the min¬ 
imal bandwidth, minimal spectral radius and minimal 
ABTE operators is listed in table lVllI 

As in the D 4 _ 3 case, the minimum spectral radius 
operator has a much smaller spectral radius than the 
other ones but, again, we did not manage to stabilize 
it with dissipation (at least, with reasonable amounts of 
it) in the 3D case with non-constant coefficients. The er¬ 
rors for the minimum ABTE operator are significantly 
smaller than the minimum bandwidth one, which is re¬ 
flected in practice by smaller errors; in addition, the op¬ 
erator could be stabilized with significantly less artificial 
dissipation. 

It should be noted at this point that we did not man¬ 
age to stabilize the Dg _5 operator with a naive adapted 
Kreiss-Oliger like dissipation prescription. We tried ap¬ 
plying Kreiss-Oliger dissipation in the interior of the 
domain, and applying no dissipation near the bound¬ 
ary where the centered stencils cannot be applied. This 
failed. Presumably this is caused by the fact that the 
dissipation operator is not negative semi-definite with 
respect to the SBP norm. Only after constructing bound¬ 
ary dissipation operators following the approach of Ref. 
118j], did we arrive at a stable scheme involving the Dg _5 
operator. 

The family of Dg _7 operators has five independent pa¬ 
rameters. When attempting to obtain minimum band¬ 
width operators by setting the coefficients 9 / 1 , 9 , 9 / 2 , 11 / 
9 / 2 , 12 / 1 / 2 , 13 / and 1 / 3,13 to zero, we numerically find 24 solu¬ 
tions, of which 16 are complex and 8 are real. However, 
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none of these solutions yields a positive definite norm. 
The minimum spectral radius operator has so large er¬ 
ror coefficients that we could not stabilize it with reason¬ 
able amounts of dissipation in the 3D non-constant co¬ 
efficient case. The minimum ABTE operator has a spec¬ 
tral radius larger than 60000 and so would require very 
small time steps when explicit time evolution is used. 

Since none of the operators considered so far was us¬ 
able in practice, we experimented with several other 
ways of choosing the parameters. First, we tried to 
minimize a weighted average of the spectral radius and 
ABTE for different weight values, but none of these op¬ 
erators turned out to be useful, even though their prop¬ 
erties were much improved. By minimizing the sum of 
the squares of the difference between error coefficients 
in neighboring boundary points, we next attempted to 
reduce the noise produced in the boundary region. Even 
when weighted with the spectral radius, these opera¬ 
tors proved not to be an improvement. Finally, based 
on the observation that some choices of parameters lead 
to very large values in the inverse of the norm, which 
directly affects the dissipation operator near the bound¬ 
ary (see Eq. 11 OH , we speculated that for some parameter 
sets the dissipation operator might make things worse 
near the boundary, and experimented with choosing pa¬ 
rameters that would minimize the condition number of 
the norm in combination with any of the previously 
mentioned properties. 

However, even though we were able to find param¬ 
eter sets that looked reasonable with respect to spec¬ 
tral radius, error coefficients and properties of the norm, 
none of the operators that we constructed could be sta¬ 
bilized with any amount of dissipation in the 3D non¬ 
constant coefficient case. This is presumably related 
to some important properties not hol ding in the non¬ 
diagonal case, as discussed in Section Til Bl Of course, 
it could still happen that usable Dg_y operators do exist, 
using some other criteria to choose parameters. 


2. Numerical tests 

For the restricted full operators we usually have to 
add dissipation. There are several causes for this, which 
are well understood; we have reviewed them in Section 

gh 

a. The operator D 4 _ 3 . One main driving point be¬ 
hind using operators based on non-diagonal norms is 
that their order of accuracy near the boundary is higher. 
Our operators based on restricted full norms drop one 
order of accuracy near the boundary, which means that 
the global convergence order is, in theory, not affected. 
(See the standard Dg _3 operator, described in Section 
IV A 2 al and illustrated in Figure|21 for an example where 
this is false in practice.) 

Figure [ 6 ] shows the results of convergence tests in the 
Loo norm for D 4 _3 operators, for both the minimum 
bandwidth and the new, optimized, minimum ABTE 


Minimum bandwidth 4-3 operators 



Time 


Minimum ABTE 4-3 operators 



Figure 6: Convergence exponents for the minimum bandwidth 
(top) and the minimum ABTE (bottom) D 4 _3 operators. 


operator. In both cases the global convergence exponent 
is very close to four. Different from the operators based 
on diagonal norms, the convergence order is also almost 
constant in time. This may be so because the bound¬ 
ary is here not the main cause of discretization error, so 
that the resulting accuracy is here independent of what 
kind of feature of the solution is currently propagating 
through the boundary. 

Figure |Z| compares the two D 4 _ 4 operators at a fixed, 
high resolution. The two operators lead to very similar 
Loo norms of the solution error. This difference in accu¬ 
racy is much smaller than it was for the operators based 
on diagonal norms. Altogether, the optimization leads 
to no practical advantage. 

b. The operator Dg_ 5 . The operator D^-s is expected 
to have the highest global order of accuracy of all the op¬ 
erators discussed in this paper (since we were not able 
to stabilize the Ds _7 operator with dissipation). We do 
in fact see sixth order convergence, but only when using 
a sufficiently accurate time integration scheme. Figure 
[8] shows the results of convergence tests in the Leo norm 
for the new, optimized D 6_5 operator with dissipation, 
using a Courant factor A = 0.25, for both fourth and 
sixth order accurate Runge-Kutta time integrators (RK4 
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4-3 operators 



Time 


Figure 7: Accuracy comparison of the two types of D 4 _ 3 oper¬ 
ators shown in the previous figure in the Lnorm. Contrary 
to the results for the Dg_ 3 and Dg 4 , here the difference in errors 
is minimal. 


and RK 6 ). 

When RK4 is used, the convergence exponent drops 
from about 6 to close to 5 for our highest resolution. 
This effect is not present when we use RK 6 , indicating 
that it is the time integrator's lower convergence order 
that actually poisons the results when resolution is in¬ 
creased. Instead of using a higher order time integrator, 
it would also have been possible to reduce the time step 
size. Especially in complicated geometries, an adaptive 
step size control is very convenient; this lets one spec¬ 
ify the desired time integration error, and the step size is 
automatically adjusted accordingly. 


C. Comparison between operators 


1. Comparison of accuracy 


We now compare the operators based on a diagonal 
and on a restricted full norm that we described and ex¬ 
amined above. Figure[9]shows, for two different resolu¬ 
tions, the solution errors for all our new operators. As 
one can see, our best performing operators are Ds _4 and 
D 6 _ 5 . One can also see that, for the highest resolution 
shown there, which corresponds to 161 3 grid points per 
block, there is a difference of five orders of magnitude be¬ 
tween the errors of Dg _5 and D 3 _i . This demonstrates 
nicely the superiority of our new high order operators 
when a high accuracy is desired. 

Even for the lowest resolution shown here, which uses 
only 21 3 grid points per block, the difference is still more 
than one order of magnitude, indicating that Dg _5 or 
D 4 _3 would be fine choices there. 


Minimum ABTE 6-5 operators (RK4) 



Time 


Minimum ABTE 6-5 operators (RK6) 



Time 


Figure 8: Convergence for the optimized Dg _5 operator with 
dissipation and Runge-Kutta time integrators of order four 
(top) and six (bottom), respectively. We see a slightly lower 
convergence order for the highest resolution when RK4 is 
used. This effect is not present with the RK6 integrator. The 
lower convergence order indicates that the accuracy of the spa¬ 
tial finite differencing operators is high enough, so that the 
overall error is dominated by the accuracy of the time inte¬ 
grator. 


2. Comparison of cost 

Higher order operators allow higher accuracy with 
the same number of grid points, but they also require 
more operations per grid point and thus have a higher 
cost. Table IVim compares actual run times per time in¬ 
tegrator step for runs with 161 3 grid points per block. 
These measured costs include not only calculating the 
right hand side, which requires taking the derivatives, 
but also applying the boundary conditions, which re¬ 
quires decomposing the system into its characteristic 
modes, and also some inter-processor communication. 
The time spent in the time integration itself is negligi¬ 
ble. 

Larger stencils increase the cost slightly. This increase 
in cost is more than made up by the increase in accu¬ 
racy, as shown above. For operators based on diago¬ 
nal norms, adding dissipation to the system increases 
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All operators (dx=0.1) 



All operators (dx=0.0125) 



Errors at fixed time t=4.4 



Figure 9: Comparison of the errors for all the new, optimized, 
differencing operators constructed in this paper. The top (mid¬ 
dle) plot shows a comparison of all the unique and optimized 
operators at low (high) resolutions. The bottom plot shows the 
error at t = 4.4 for all resolutions. The most successful opera¬ 
tors are the optimized D 5 _ 5 , D 4 _ 3 , and D§_ 4 . 


the cost only marginally. The dissipation operators for 
derivatives based on non-diagonal norms are more com¬ 
plex to calculate and increase the run time noticeably. 
All numbers were obtained from a straightforward legi¬ 
ble implementation in Fortran, without spending much 
effort on optimizing the code for performance. 

The effect of higher order operators on the overall run 


Table VIII: Relative operator costs for our test problem, mea¬ 
sured in seconds per time step with an RK4 integrator. The 
number in parentheses for the Dg _5 used an RK6 time integra¬ 
tor, which requires seven iterations instead of four. The num¬ 
bers marked with an asterisk were obtained from runs using a 
larger number of processors and are probably larger than they 
would be otherwise due to scaling inefficiencies. 


Operator 

Time [sec] 
without dissipation 

Time [sec] 
with dissipation 

D 2 _i 

26.9 

27.7 

D 4 _2 

27.9 

29.7 

D 6-3 

30.5 

32.0 

£>8-4 

34.4 

35.3 

D 4-3 

28.6 

44.8 

£*6-5 

41.6* 

39.1 (86.2*) 


time is not very pronounced. As they do not increase 
the storage requirements either, the only reason speak¬ 
ing against using them (for smooth problems) seems to 
be the effort one has to spend constructing and imple¬ 
menting them. 


VI. CONCLUSION 

Let us summarize the main points of this paper. 
We have explicitly constructed accurate, high-order fi¬ 
nite differencing operators which satisfy summation by 
parts. This construction is not unique, and it is necessary 
to specify some free parameters. We have considered 
several optimization criteria to define these parameters; 
namely, (a) a minimum bandwidth, (b) a minimum of 
the truncation error on the boundary points, (c) a mini¬ 
mum spectral radius, and a combination of these. 

We examined in detail a set of operators that are 
up to tenth order accurate. We found that minimum 
bandwidth operators may have a large spectral radius 
or truncation error near the boundary. Optimizing for 
these two criteria can, surprisingly, reduce the opera¬ 
tors' spectral radius by three orders and their accuracy 
near the boundary by two orders of magnitude. 

Some of the finite differencing operators require ar¬ 
tificial dissipation. We have therefore also constructed 
high-order dissipation operators, compatible with the 
above finite differencing operators, and also semi- 
definite with respect to the SBP scalar product. 

We tested the stability and accuracy of these opera¬ 
tors by evolving a scalar wave equation on a spherical 
domain. Our domain is split into seven blocks, each 
discretized with a structured grid. These blocks are 
connected through penalty boundary conditions. We 
demonstrated that the optimized finite differencing op¬ 
erators have also far superior properties in practice. The 
most accurate operators are Dg _4 and Dg_ 5 ; the latter 
is in our setup five orders of magnitude more accurate 
than a simple D 2 _i operator. 
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Appendix A: OPERATOR COEFFICIENTS 

We provide, for the reader's convenience, the coeffi¬ 
cients for the derivative and dissipation operators that 
we constructed above. Since the values of these coeffi¬ 
cients themselves are only of limited interest, and since 
there is a large danger of introducing errors when type¬ 
setting these coefficients, we make them available elec¬ 
tronically instead. We also make a Cactus M thorn 
SummationByParts available via anonymous CVS. This 
thorn implements the derivative and dissipation oper¬ 
ators. Our web pages |^9] contain instructions for ac¬ 
cessing these. For the sake of continuity, we will also 
make the coefficients available on www.arxiv.org to¬ 
gether with this article. 

We distribute the coefficients as a set of files, where 
each file defines on operator. The content of the file is 
written in a Fortran-like pseudo language that defines 
the coefficients in declarations like 

a(l) = 0.5 


and sometimes makes use of additional constants, as in 
xl = 3 

a(2) = xl + 4 

We write here a(l) as ci\ and q(2,3) as < 723 . 


1. Derivative operators 

The derivative operators D 2 -I, D 4 - 2 , D 6 _ 3 , Dg_ 4 , 
D 4 _ 3 , and D 6 _ 5 are defined via coefficients fl, and q r j. 
In the interior of the domain it is 

1 s 

DijUj = — E a j ( u i+j — u i—j) (Al) 

7=1 

and near the left boundary it is (i.e. i = 1 , b ) 

1 s 

DijUj = — Y_J (jjiUj- (A2) 

j=1 

At the right boundary the same coefficients are used in 
opposite order and with opposite sign. 


2. Dissipation operators 

a. Dissipation operators based on diagonal norms 


The dissipation operators corresponding to D 2 _i, 
D 4 _ 2 , Dg_ 3 , and Dg _4 are defined via coefficients a,j and 
Hi- 

In the interior of the domain it is 


A-ijUj 


e 


qoUi + E <7; («/-;• + Ui+j) 
7=1 


and near the boundary it is 


AijUj 


£ 

2 


E a ji U j' 

7=1 


(A3) 


(A4) 


where e > 0 selects the amount of dissipation and is 
usually of order unity. 


b. Dissipation operators based on non-diagonal norms 

The dissipation operators corresponding to 04.3 and 
Dg_ 5 ... are more complicated, since they depend on the 
user parameters specifying the number of grid points, 
N, and the size of the transition region (i.e. the region 
where Bp is different from 1 ). 
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The dissipation operators are then constructed ac¬ 
cording to Eq. flOl 

A 2p = -^ h2p Z~ lD J,B p D p . (A5) 

The coefficients for the inverse of the norm, E -1 , are 
provided in the boundary region only. In the files this 
inverse is denoted by sigmafi, j). In the interior the 
norm (and its inverse) is diagonal with value 1 . 

In the D 4_3 case, the N x N matrix D 2 defining the 
consistent approximation of d 2 /dx 2 is given by 


/I -2 10 

1-2 10 
0 1-2 1 



V 


\ 


1-2 10 
0 1-2 1 

0 1 - 21 / 


(A 6 ) 


while the diagonal matrix £>2 has the value h at either 
boundary and increases linearly to the value 1 across the 
user defined transition region. 

In the Dj ,_5 case the matrix defining the consistent ap¬ 
proximation of d?/dy? near the left boundary D 3 is 


D 


l 

3 



3 

3 

3 

1 


-3 1 0 

-3 1 0 

-3 1 0 

3 -3 1 


V 


\ 

J- 


(A 7) 


while at the right boundary K is 


D5 


/ 

1 

¥ 

V 


-1 3-3 10 

0-1 3-31 

0-1 3-31 

0-1 3-31/ 


(A 8 ) 


Since the values on the diagonal in the interior of D( 
and Dj are —3 and 3, respectively, it is impossible to 
construct a single matrix D 3 to cover the whole do¬ 
main. However, since both matrix products {D^) T D{ ; 
and (D' 3 ) t D' 3 result in the same interior operator, dissi¬ 
pation operators can be constructed in the left and right 
domain separately and then combined into a global op¬ 
erator. The diagonal matrix B 3 has the values h 2 at 
the boundary and and 1 in the interior, and a third or¬ 
der polynomial with zero derivative at either end of the 
transition region is used to smoothly connect the bound¬ 
ary with the interior. 
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